github link: https://github.com/becca-belmonte/H99_RNAseq
First, we can see that the samples primarily differ in expression due to being control or H99 samples. We can also see the top differentially expressed genes, both those higher in control and those higher in H99 mutants.
#rownames(col_data) %in% colnames(raw_counts)
matrix_raw_counts <- as.matrix(raw_counts)
dds <- DESeqDataSetFromMatrix(countData = matrix_raw_counts,
colData = col_data,
design = ~ Group)
dds <- estimateSizeFactors(dds)
vsdB = varianceStabilizingTransformation(dds)
plotPCA(vsdB,intgroup = c("Group"))
ddsTC <- DESeq(dds)
resTC <- results(ddsTC)
resTC$gene <- rownames(resTC)
tab_resTC <- as.data.frame(resTC)
tab_resTC$gene <- rownames(tab_resTC)
tab_resTC <- cbind(tab_resTC, Gene_ID_Flybase[match(tab_resTC$gene, Gene_ID_Flybase$FBgn_ID), "Gene_name"])
colnames(tab_resTC)[ncol(tab_resTC)] = "Gene_name"
tab_resTC <- cbind(tab_resTC, Gene_ID_Flybase[match(tab_resTC$gene, Gene_ID_Flybase$FBgn_ID), "CG_ID"])
colnames(tab_resTC)[ncol(tab_resTC)] = "CG_ID"
tab_resTC <- tab_resTC %>%
mutate(Gene_name = coalesce(Gene_name, gene))
tab_top <- tab_resTC %>%
filter(abs(log2FoldChange) > 1 & padj < 0.05) %>%
arrange(padj)
tab_top_down <- tab_top %>%
filter(log2FoldChange < 0)
kable(head(tab_top, 10)) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
scroll_box(width = "100%")
| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | gene | Gene_name | CG_ID | |
|---|---|---|---|---|---|---|---|---|---|
| FBgn0085466 | 43263.1941 | 4.824477 | 0.2282038 | 21.14109 | 0 | 0 | FBgn0085466 | FBgn0085466 | NA |
| FBgn0020493 | 4746.8218 | -5.549997 | 0.2652973 | -20.91992 | 0 | 0 | FBgn0020493 | Dad | CG5201 |
| FBgn0037930 | 6747.1360 | 3.919628 | 0.2287028 | 17.13852 | 0 | 0 | FBgn0037930 | CG14715 | CG14715 |
| FBgn0052085 | 1674.8032 | 8.845837 | 0.5171809 | 17.10395 | 0 | 0 | FBgn0052085 | FBgn0052085 | NA |
| FBgn0023129 | 25753.8140 | 2.400866 | 0.1460641 | 16.43707 | 0 | 0 | FBgn0023129 | aay | CG3705 |
| FBgn0025678 | 53362.0963 | -3.282450 | 0.2067646 | -15.87530 | 0 | 0 | FBgn0025678 | CaBP1 | CG5809 |
| FBgn0037222 | 22074.0559 | -2.704999 | 0.1858704 | -14.55314 | 0 | 0 | FBgn0037222 | CG14642 | CG14642 |
| FBgn0039099 | 15792.9504 | -2.203793 | 0.1717471 | -12.83161 | 0 | 0 | FBgn0039099 | GILT2 | CG10157 |
| FBgn0051813 | 3122.6405 | 3.710390 | 0.2991312 | 12.40389 | 0 | 0 | FBgn0051813 | nur | CG31813 |
| FBgn0267472 | 396.0726 | -5.570781 | 0.4826268 | -11.54263 | 0 | 0 | FBgn0267472 | asRNA:CR45822 | CR45822 |
#write.csv(tab_top, file = "Results/tab_top.csv")
To get a better idea of what these genes are doing, I did a gene set enrichment analysis. First I looked at the reactome and unfortunately, none of them are significant with the adjusted p-value. But, we do see that Golgi Associated Vesicle Biogenesis is particularly enriched in H99 mutants when compared to the control. I then looked at KEGG pathways and saw that the phagosome pathway was enriched in control samples when compared to H99 mutants.
colors_vol <- c(rev(viridis(2)), "darkgrey")
ranks <- tab_resTC$log2FoldChange[!is.na(tab_resTC$log2FoldChange)]
names(ranks) <- tab_resTC$CG_ID[!is.na(tab_resTC$log2FoldChange)]
#head(ranks)
#barplot(sort(ranks, decreasing = T))
fgseaRes <- fgsea(t_Drosophila_REACTOME, ranks, minSize = 15, maxSize = 500)
head(fgseaRes[order(padj, -abs(NES)), ], n=10)
## pathway
## 1: Golgi Associated Vesicle Biogenesis
## 2: Phosphorylation of PER and TIM
## 3: Signaling by Non-Receptor Tyrosine Kinases
## 4: Signaling by PTK6
## 5: Intra-Golgi traffic
## 6: RNA Polymerase II Promoter Escape
## 7: RNA Polymerase II Transcription Initiation
## 8: RNA Polymerase II Transcription Initiation And Promoter Clearance
## 9: RNA Polymerase II Transcription Pre-Initiation And Promoter Opening
## 10: MAP2K and MAPK activation
## pval padj log2err ES NES size
## 1: 0.002945886 1 0.4317077 0.7107300 1.724223 27
## 2: 0.020256094 1 0.3524879 0.7349901 1.554007 15
## 3: 0.028444229 1 0.3524879 0.6875806 1.527877 18
## 4: 0.028444229 1 0.3524879 0.6875806 1.527877 18
## 5: 0.022954825 1 0.3524879 0.6294295 1.526989 27
## 6: 0.021969374 1 0.3524879 -0.6054070 -1.508175 38
## 7: 0.021969374 1 0.3524879 -0.6054070 -1.508175 38
## 8: 0.021969374 1 0.3524879 -0.6054070 -1.508175 38
## 9: 0.021969374 1 0.3524879 -0.6054070 -1.508175 38
## 10: 0.070981211 1 0.2450418 0.6082707 1.448003 24
## leadingEdge
## 1: CG32683,CG4349,CG5711
## 2: CG3772,CG2048
## 3: CG3927,CG10384,CG3875,CG4021
## 4: CG3927,CG10384,CG3875,CG4021
## 5: CG10592,CG5150,CG3292,CG7980,CG8105,CG8147
## 6: CG14718,CG7562,CG9879,CG3180
## 7: CG14718,CG7562,CG9879,CG3180
## 8: CG14718,CG7562,CG9879,CG3180
## 9: CG14718,CG7562,CG9879,CG3180
## 10: CG32683,CG5711,CG31832,CG7524,CG8642,CG10359,...
plotEnrichment(t_Drosophila_REACTOME[["Golgi Associated Vesicle Biogenesis"]], ranks) +
ggtitle("Golgi Associated Vesicle Biogenesis")
fgseaRes <- fgseaRes %>%
mutate(color = factor(case_when(NES >= 1 & pval < 0.05 ~ "blue",
NES <= -1 & pval < 0.05 ~ "red",
NES %in% c(-1:1) ~ "grey")))
ggplot(fgseaRes, aes(x = NES, y = -log10(pval))) +
geom_point(aes(color = color))+
geom_text_repel(data = (fgseaRes %>% filter(pval < 0.025)),aes(label = pathway)) +
ggtitle("Enriched pathways without regulators of apoptosis") +
theme_light()+
scale_colour_manual(limits=c("blue","red", "black"),
values = colors_vol,
labels=c("Enriched in H99","Enriched in control","ns")) +
labs(color = "")
fgseaRes_KEGG <- fgsea(t_Drosophila_KEGG, ranks, minSize = 15, maxSize = 500)
head(fgseaRes_KEGG[order(padj, -abs(NES)), ], n=10)
## pathway pval padj log2err
## 1: Phagosome 0.04447955 1 0.3217759
## 2: Biosynthesis of unsaturated fatty acids 0.10052910 1 0.1864326
## 3: Galactose metabolism 0.14010508 1 0.1552420
## 4: Basal transcription factors 0.11981567 1 0.1957890
## 5: Insect hormone biosynthesis 0.14712644 1 0.1752040
## 6: Apoptosis_multiple species 0.17040359 1 0.1596467
## 7: ABC transporters 0.19964974 1 0.1275053
## 8: Cysteine and methionine metabolism 0.19237435 1 0.1294429
## 9: N Glycan biosynthesis 0.17667845 1 0.1372508
## 10: Longevity regulating pathway 0.16067146 1 0.1709323
## ES NES size
## 1: -0.5395880 -1.446518 57
## 2: -0.6490166 -1.365043 16
## 3: -0.6023040 -1.314794 19
## 4: 0.5035614 1.293475 32
## 5: 0.5763314 1.264763 16
## 6: 0.5861498 1.264691 15
## 7: -0.5730456 -1.250924 19
## 8: -0.5299475 -1.243260 27
## 9: -0.5038413 -1.228812 34
## 10: 0.4696273 1.223888 38
## leadingEdge
## 1: CG8310,CG1076,CG9906,CG1924,CG7794,CG7678,...
## 2: CG15531,CG9747,CG18609,CG10849
## 3: CG8695,CG14934,CG12766,CG32444,CG10638,CG5165,...
## 4: CG15632,CG11639,CG5444,CG1276
## 5: CG13478,CG40486,CG31075,CG40485,CG10594,CG6578,...
## 6: CG8238,CG6531,CG8091,CG7788,CG6829,CG14902
## 7: CG4562,CG8799,CG8908,CG10181,CG31792
## 8: CG11899,CG6287,CG31115
## 9: CG4871,CG17173
## 10: CG14173,CG14167,CG7756,CG7978,CG5948,CG1506,...
plotEnrichment(t_Drosophila_KEGG[["Phagosome"]], ranks) +
ggtitle("Phagosome")
fgseaRes_KEGG <- fgseaRes_KEGG %>%
mutate(color = factor(case_when(NES >= 1 & pval < 0.05 ~ "blue",
NES <= -1 & pval < 0.05 ~ "red",
NES %in% c(-1:1) ~ "grey")))
ggplot(fgseaRes_KEGG, aes(x = NES, y = -log10(pval))) +
geom_point(aes(color = color))+
geom_text_repel(data = (fgseaRes_KEGG %>% filter(pval < 0.1)),aes(label = pathway)) +
ggtitle("Enriched KEGG pathways without regulators of apoptosis") +
theme_light()+
scale_colour_manual(limits=c("blue","red", "black"),
values = colors_vol,
labels=c("Enriched in H99","Enriched in control","ns")) +
labs(color = "")
phago_genes <- tab_resTC %>%
filter(CG_ID %in% t_Drosophila_KEGG[["Phagosome"]])
top_phago_genes <- c("FBgn0264077", "FBgn0003884", "FBgn0005671", "FBgn0040377", "FBgn0261797")
I wanted to get a better idea of which genes were lost in H99 mutants, and particularly those in the phagosome pathway.
best12 <- head(tab_top_down,12)
best12_genes <- best12$gene
col_data <- col_data %>%
arrange(ID)
normalized <- normalized_counts %>%
mutate(gene = rownames(normalized_counts)) %>%
filter(gene %in% best12_genes) %>%
select(-gene)
normalized <- as.data.frame(t(normalized))
normalized <- normalized %>%
mutate(group = rownames(normalized)) %>%
arrange(group) %>%
cbind(col_data) %>%
pivot_longer(best12_genes, names_to = "gene", values_to = "counts")
normalized$Rep <- as.factor(normalized$Rep)
normalized <- cbind(normalized, Gene_ID_Flybase[match(normalized$gene, Gene_ID_Flybase$FBgn_ID), "Gene_name"])
colnames(normalized)[ncol(normalized)] = "Gene_name"
normalized <- normalized %>%
mutate(Gene_name = coalesce(Gene_name, gene))
(rawcounts <- ggplot(normalized, aes(x = Group, y = counts)) +
geom_point(aes(color = Group), position = position_dodge(width = 0.75)) +
scale_color_manual(values = (viridis(2))) +
facet_wrap(~Gene_name, scale = "free_y") +
theme_classic() +
xlab("Genotype") +
ylab("Normalized counts") +
ggtitle("Top differentially expressed genes without apoptosis"))
col_data <- col_data %>%
arrange(ID)
normalized <- normalized_counts %>%
mutate(gene = rownames(normalized_counts)) %>%
filter(gene %in% top_phago_genes) %>%
select(-gene)
normalized <- as.data.frame(t(normalized))
normalized <- normalized %>%
mutate(group = rownames(normalized)) %>%
arrange(group) %>%
cbind(col_data) %>%
pivot_longer(top_phago_genes, names_to = "gene", values_to = "counts")
normalized$Rep <- as.factor(normalized$Rep)
normalized <- cbind(normalized, Gene_ID_Flybase[match(normalized$gene, Gene_ID_Flybase$FBgn_ID), "Gene_name"])
colnames(normalized)[ncol(normalized)] = "Gene_name"
normalized <- normalized %>%
mutate(Gene_name = coalesce(Gene_name, gene))
(rawcounts <- ggplot(normalized, aes(x = Group, y = counts)) +
geom_point(aes(color = Group), position = position_dodge(width = 0.75)) +
scale_color_manual(values = (viridis(2))) +
facet_wrap(~Gene_name, scale = "free_y") +
theme_classic() +
xlab("Genotype") +
ylab("Normalized counts") +
ggtitle("Top differentially expressed phagosome-associated genes without apoptosis"))
With the volcano plot, we can see that Dad and CaBP1 are downregulated in H99 mutants, while FBgn0085466 (which is CG34437) and CG14715 are upregulated in H99 mutants.
colors_vol <- c(rev(viridis(2)), "darkgrey")
tab_resTC <- tab_resTC %>%
mutate(color = factor(case_when(log2FoldChange >= 2 & padj < 0.05 ~ "up",
log2FoldChange <= -2 & padj < 0.05 ~ "down",
log2FoldChange %in% c(-2:2) ~ "ns"))) %>%
mutate(significance = case_when(log2FoldChange >= 2 ~ "yes", log2FoldChange <= -2 ~ "yes"))
significant <- subset(tab_resTC, significance=="yes")
sig_genes <- significant$gene
tab_resTC$color <- as.character(tab_resTC$color)
tab_resTC$color <- replace_na(tab_resTC$color,"ns")
tab_resTC <- tab_resTC %>%
mutate(Gene_name = coalesce(Gene_name, gene))
(carcass <- ggplot(tab_resTC, aes(x = log2FoldChange, y = -log10(pvalue), label = Gene_name)) +
geom_point(aes(color = color), alpha = 0.5)+
geom_text_repel(data = tab_resTC %>% filter(log2FoldChange < -2 | log2FoldChange > 2.5) %>% filter(pvalue <10^-5), aes(label = Gene_name)) +
geom_vline(xintercept=2,color="darkgrey", linetype = "dotted")+
geom_vline(xintercept=-2,color="darkgrey", linetype = "dotted")+
geom_hline(yintercept=-log10(0.05),color="darkgrey",linetype="dashed") +
theme_bw() +
scale_colour_manual(limits=c("up","down", "ns"),
values = colors_vol,
labels=c("Upregulated in H99","Upregulated in control","ns")) +
labs(color = "") +
ggtitle("Differential expression in control and H99 mutants"))
ggplotly(carcass)
The first table shows GO terms that are enriched in upregulated genes, so those that are higher in H99 mutants than control. The second table shows GO terms enriched in downregulated genes, those that are lower in H99 mutants.
mart <- useDataset(dataset = "dmelanogaster_gene_ensembl",
mart = useMart("ENSEMBL_MART_ENSEMBL"))
resultTable <- getBM(attributes = c("flybase_gene_id", "external_gene_name","go_id", "name_1006", "definition_1006"),
mart = mart)
resultTable <- resultTable[resultTable$go_id != '',]
geneID2GO <- by(resultTable$go_id,
resultTable$external_gene_name,
function(x) as.character(x))
datRef <- tab_resTC$Gene_name #full set of genes in your analysis (ie. rownames of your summarised experiment input into DEseq2)
tab_top_up <- tab_top %>%
filter(log2FoldChange > 0)
theGenes <- tab_top_up$Gene_name#set of flybase IDs that you want to query
geneNames <- tab_resTC$Gene_name
myInterestingGenes <- tab_top_up$Gene_name
geneList <- tab_resTC %>%
filter(Gene_name %in% myInterestingGenes) %>%
na.omit(geneList)
geneList <- geneList %>%
arrange(pvalue)
geneList_pvalue <- geneList[,5]
geneList_name <- as.character(geneList$Gene_name)
names(geneList_pvalue) <- geneList_name
geneList <- geneList_pvalue
#write.csv(geneList_F_UC_8,"geneList_F_UC_8.csv",col.names=F,row.names=T)
all_genes <- sort(unique(as.character(resultTable$external_gene_name)))
int_genes <- factor(as.integer(all_genes %in% geneList_name))
names(int_genes) = all_genes
GOdata <- new("topGOdata", ontology = "BP", allGenes = int_genes,
annot = annFUN.gene2GO, gene2GO = geneID2GO)
GO_results_classic <- runTest(GOdata, algorithm = "classic", statistic = "Fisher")
GO_results_elim <- runTest(GOdata, algorithm = "elim", statistic = "Fisher")
GO_results_tab_up <- GenTable(object = GOdata, classic = GO_results_classic, elim = GO_results_elim, orderBy = "elim", ranksOf = "classic", topNodes = 50)
kable(GO_results_tab_up) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>% scroll_box(width = "100%")
| GO.ID | Term | Annotated | Significant | Expected | Rank in classic | classic | elim |
|---|---|---|---|---|---|---|---|
| GO:0098754 | detoxification | 87 | 4 | 0.55 | 1 | 0.0023 | 0.0023 |
| GO:0043153 | entrainment of circadian clock by photop… | 12 | 2 | 0.08 | 2 | 0.0025 | 0.0025 |
| GO:0009785 | blue light signaling pathway | 1 | 1 | 0.01 | 4 | 0.0064 | 0.0064 |
| GO:0010835 | regulation of protein ADP-ribosylation | 1 | 1 | 0.01 | 5 | 0.0064 | 0.0064 |
| GO:0014853 | regulation of excitatory postsynaptic me… | 1 | 1 | 0.01 | 6 | 0.0064 | 0.0064 |
| GO:0018312 | peptidyl-serine ADP-ribosylation | 1 | 1 | 0.01 | 7 | 0.0064 | 0.0064 |
| GO:0030216 | keratinocyte differentiation | 1 | 1 | 0.01 | 8 | 0.0064 | 0.0064 |
| GO:0032481 | positive regulation of type I interferon… | 1 | 1 | 0.01 | 9 | 0.0064 | 0.0064 |
| GO:0050843 | S-adenosylmethionine catabolic process | 1 | 1 | 0.01 | 10 | 0.0064 | 0.0064 |
| GO:0050980 | detection of light stimulus involved in … | 1 | 1 | 0.01 | 11 | 0.0064 | 0.0064 |
| GO:0071465 | cellular response to desiccation | 1 | 1 | 0.01 | 12 | 0.0064 | 0.0064 |
| GO:1901054 | sarcosine biosynthetic process | 1 | 1 | 0.01 | 13 | 0.0064 | 0.0064 |
| GO:0030001 | metal ion transport | 285 | 6 | 1.82 | 25 | 0.0096 | 0.0096 |
| GO:0007610 | behavior | 588 | 9 | 3.75 | 26 | 0.0122 | 0.0122 |
| GO:0006111 | regulation of gluconeogenesis | 2 | 1 | 0.01 | 27 | 0.0127 | 0.0127 |
| GO:0009588 | UV-A, blue light phototransduction | 2 | 1 | 0.01 | 28 | 0.0127 | 0.0127 |
| GO:0009639 | response to red or far red light | 2 | 1 | 0.01 | 29 | 0.0127 | 0.0127 |
| GO:0010114 | response to red light | 2 | 1 | 0.01 | 30 | 0.0127 | 0.0127 |
| GO:0015770 | sucrose transport | 2 | 1 | 0.01 | 31 | 0.0127 | 0.0127 |
| GO:0048150 | behavioral response to ether | 2 | 1 | 0.01 | 32 | 0.0127 | 0.0127 |
| GO:0061965 | positive regulation of entry into reprod… | 2 | 1 | 0.01 | 33 | 0.0127 | 0.0127 |
| GO:0071000 | response to magnetism | 2 | 1 | 0.01 | 34 | 0.0127 | 0.0127 |
| GO:0071929 | alpha-tubulin acetylation | 2 | 1 | 0.01 | 35 | 0.0127 | 0.0127 |
| GO:0120176 | positive regulation of torso signaling p… | 2 | 1 | 0.01 | 36 | 0.0127 | 0.0127 |
| GO:0001752 | compound eye photoreceptor fate commitme… | 77 | 3 | 0.49 | 39 | 0.0130 | 0.0130 |
| GO:0042706 | eye photoreceptor cell fate commitment | 77 | 3 | 0.49 | 40 | 0.0130 | 0.0130 |
| GO:0007218 | neuropeptide signaling pathway | 78 | 3 | 0.50 | 41 | 0.0135 | 0.0135 |
| GO:0048878 | chemical homeostasis | 314 | 6 | 2.00 | 42 | 0.0149 | 0.0149 |
| GO:0042592 | homeostatic process | 510 | 8 | 3.25 | 43 | 0.0156 | 0.0156 |
| GO:0046552 | photoreceptor cell fate commitment | 84 | 3 | 0.54 | 45 | 0.0164 | 0.0164 |
| GO:0098659 | inorganic cation import across plasma me… | 32 | 2 | 0.20 | 46 | 0.0176 | 0.0176 |
| GO:0099587 | inorganic ion import across plasma membr… | 32 | 2 | 0.20 | 47 | 0.0176 | 0.0176 |
| GO:0030534 | adult behavior | 157 | 4 | 1.00 | 48 | 0.0179 | 0.0179 |
| GO:0006564 | L-serine biosynthetic process | 3 | 1 | 0.02 | 50 | 0.0190 | 0.0190 |
| GO:0008302 | female germline ring canal formation, ac… | 3 | 1 | 0.02 | 51 | 0.0190 | 0.0190 |
| GO:0030046 | parallel actin filament bundle assembly | 3 | 1 | 0.02 | 52 | 0.0190 | 0.0190 |
| GO:0045472 | response to ether | 3 | 1 | 0.02 | 53 | 0.0190 | 0.0190 |
| GO:0046498 | S-adenosylhomocysteine metabolic process | 3 | 1 | 0.02 | 54 | 0.0190 | 0.0190 |
| GO:0007155 | cell adhesion | 242 | 5 | 1.54 | 57 | 0.0190 | 0.0190 |
| GO:0008543 | fibroblast growth factor receptor signal… | 37 | 2 | 0.24 | 59 | 0.0232 | 0.0232 |
| GO:0043270 | positive regulation of ion transport | 37 | 2 | 0.24 | 60 | 0.0232 | 0.0232 |
| GO:0044344 | cellular response to fibroblast growth f… | 37 | 2 | 0.24 | 61 | 0.0232 | 0.0232 |
| GO:0071774 | response to fibroblast growth factor | 37 | 2 | 0.24 | 62 | 0.0232 | 0.0232 |
| GO:0008340 | determination of adult lifespan | 172 | 4 | 1.10 | 63 | 0.0241 | 0.0241 |
| GO:0000768 | syncytium formation by plasma membrane f… | 38 | 2 | 0.24 | 64 | 0.0243 | 0.0243 |
| GO:0006949 | syncytium formation | 38 | 2 | 0.24 | 65 | 0.0243 | 0.0243 |
| GO:0007520 | myoblast fusion | 38 | 2 | 0.24 | 66 | 0.0243 | 0.0243 |
| GO:0014902 | myotube differentiation | 38 | 2 | 0.24 | 67 | 0.0243 | 0.0243 |
| GO:0140253 | cell-cell fusion | 38 | 2 | 0.24 | 68 | 0.0243 | 0.0243 |
| GO:0048663 | neuron fate commitment | 99 | 3 | 0.63 | 70 | 0.0253 | 0.0253 |
tab_top_down <- tab_top %>%
filter(log2FoldChange < 0)
theGenes <- tab_top_down$Gene_name#set of flybase IDs that you want to query
geneNames <- tab_resTC$Gene_name
myInterestingGenes <- tab_top_down$Gene_name
geneList <- tab_resTC %>%
filter(Gene_name %in% myInterestingGenes) %>%
na.omit(geneList)
geneList <- geneList %>%
arrange(pvalue)
geneList_pvalue <- geneList[,5]
geneList_name <- as.character(geneList$Gene_name)
names(geneList_pvalue) <- geneList_name
geneList <- geneList_pvalue
#write.csv(geneList_F_UC_8,"geneList_F_UC_8.csv",col.names=F,row.names=T)
all_genes <- sort(unique(as.character(resultTable$external_gene_name)))
int_genes <- factor(as.integer(all_genes %in% geneList_name))
names(int_genes) = all_genes
GOdata <- new("topGOdata", ontology = "BP", allGenes = int_genes,
annot = annFUN.gene2GO, gene2GO = geneID2GO)
GO_results_classic <- runTest(GOdata, algorithm = "classic", statistic = "Fisher")
GO_results_elim <- runTest(GOdata, algorithm = "elim", statistic = "Fisher")
GO_results_tab_down <- GenTable(object = GOdata, classic = GO_results_classic, elim = GO_results_elim, orderBy = "elim", ranksOf = "classic", topNodes = 50)
kable(GO_results_tab_down) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>% scroll_box(width = "100%")
| GO.ID | Term | Annotated | Significant | Expected | Rank in classic | classic | elim |
|---|---|---|---|---|---|---|---|
| GO:0007094 | mitotic spindle assembly checkpoint sign… | 15 | 3 | 0.10 | 1 | 0.00013 | 0.00013 |
| GO:0007443 | Malpighian tubule morphogenesis | 44 | 4 | 0.30 | 13 | 0.00022 | 0.00022 |
| GO:0035002 | liquid clearance, open tracheal system | 24 | 3 | 0.16 | 19 | 0.00055 | 0.00055 |
| GO:0048133 | male germ-line stem cell asymmetric divi… | 8 | 2 | 0.05 | 27 | 0.00125 | 0.00125 |
| GO:0008362 | chitin-based embryonic cuticle biosynthe… | 17 | 2 | 0.12 | 44 | 0.00582 | 0.00582 |
| GO:0010991 | negative regulation of SMAD protein comp… | 1 | 1 | 0.01 | 46 | 0.00680 | 0.00680 |
| GO:0032499 | detection of peptidoglycan | 1 | 1 | 0.01 | 47 | 0.00680 | 0.00680 |
| GO:0033316 | meiotic spindle assembly checkpoint sign… | 1 | 1 | 0.01 | 48 | 0.00680 | 0.00680 |
| GO:0042940 | D-amino acid transport | 1 | 1 | 0.01 | 49 | 0.00680 | 0.00680 |
| GO:0043060 | meiotic metaphase I plate congression | 1 | 1 | 0.01 | 50 | 0.00680 | 0.00680 |
| GO:0060394 | negative regulation of pathway-restricte… | 1 | 1 | 0.01 | 51 | 0.00680 | 0.00680 |
| GO:0061734 | parkin-mediated stimulation of mitophagy… | 1 | 1 | 0.01 | 52 | 0.00680 | 0.00680 |
| GO:1901526 | positive regulation of mitophagy | 1 | 1 | 0.01 | 53 | 0.00680 | 0.00680 |
| GO:1901693 | negative regulation of compound eye reti… | 1 | 1 | 0.01 | 54 | 0.00680 | 0.00680 |
| GO:1904068 | G protein-coupled receptor signaling pat… | 1 | 1 | 0.01 | 55 | 0.00680 | 0.00680 |
| GO:1905342 | positive regulation of protein localizat… | 1 | 1 | 0.01 | 56 | 0.00680 | 0.00680 |
| GO:0006030 | chitin metabolic process | 21 | 2 | 0.14 | 68 | 0.00883 | 0.00883 |
| GO:0035214 | eye-antennal disc development | 64 | 3 | 0.44 | 69 | 0.00938 | 0.00938 |
| GO:0060438 | trachea development | 23 | 2 | 0.16 | 71 | 0.01054 | 0.01054 |
| GO:0051240 | positive regulation of multicellular org… | 197 | 5 | 1.34 | 73 | 0.01094 | 0.01094 |
| GO:0031644 | regulation of nervous system process | 24 | 2 | 0.16 | 74 | 0.01145 | 0.01145 |
| GO:0044057 | regulation of system process | 69 | 3 | 0.47 | 76 | 0.01152 | 0.01152 |
| GO:0042753 | positive regulation of circadian rhythm | 25 | 2 | 0.17 | 77 | 0.01239 | 0.01239 |
| GO:0061983 | meiosis II cell cycle process | 26 | 2 | 0.18 | 79 | 0.01337 | 0.01337 |
| GO:0006037 | cell wall chitin metabolic process | 2 | 1 | 0.01 | 81 | 0.01356 | 0.01356 |
| GO:0006038 | cell wall chitin biosynthetic process | 2 | 1 | 0.01 | 82 | 0.01356 | 0.01356 |
| GO:0007066 | female meiosis sister chromatid cohesion | 2 | 1 | 0.01 | 83 | 0.01356 | 0.01356 |
| GO:0009305 | protein biotinylation | 2 | 1 | 0.01 | 84 | 0.01356 | 0.01356 |
| GO:0010460 | positive regulation of heart rate | 2 | 1 | 0.01 | 85 | 0.01356 | 0.01356 |
| GO:0018054 | peptidyl-lysine biotinylation | 2 | 1 | 0.01 | 86 | 0.01356 | 0.01356 |
| GO:0045760 | positive regulation of action potential | 2 | 1 | 0.01 | 87 | 0.01356 | 0.01356 |
| GO:0048075 | positive regulation of eye pigmentation | 2 | 1 | 0.01 | 88 | 0.01356 | 0.01356 |
| GO:0048078 | positive regulation of compound eye pigm… | 2 | 1 | 0.01 | 89 | 0.01356 | 0.01356 |
| GO:0051971 | positive regulation of transmission of n… | 2 | 1 | 0.01 | 90 | 0.01356 | 0.01356 |
| GO:0060279 | positive regulation of ovulation | 2 | 1 | 0.01 | 91 | 0.01356 | 0.01356 |
| GO:0071110 | histone biotinylation | 2 | 1 | 0.01 | 92 | 0.01356 | 0.01356 |
| GO:0110123 | regulation of myotube cell migration | 2 | 1 | 0.01 | 93 | 0.01356 | 0.01356 |
| GO:0110125 | negative regulation of myotube cell migr… | 2 | 1 | 0.01 | 94 | 0.01356 | 0.01356 |
| GO:1904061 | positive regulation of locomotor rhythm | 2 | 1 | 0.01 | 95 | 0.01356 | 0.01356 |
| GO:1904457 | positive regulation of neuronal action p… | 2 | 1 | 0.01 | 96 | 0.01356 | 0.01356 |
| GO:1905050 | positive regulation of metallopeptidase … | 2 | 1 | 0.01 | 97 | 0.01356 | 0.01356 |
| GO:0007448 | anterior/posterior pattern specification… | 27 | 2 | 0.18 | 112 | 0.01437 | 0.01437 |
| GO:0015807 | L-amino acid transport | 28 | 2 | 0.19 | 113 | 0.01541 | 0.01541 |
| GO:1902475 | L-alpha-amino acid transmembrane transpo… | 28 | 2 | 0.19 | 114 | 0.01541 | 0.01541 |
| GO:0043171 | peptide catabolic process | 29 | 2 | 0.20 | 115 | 0.01648 | 0.01648 |
| GO:0000705 | achiasmate meiosis I | 3 | 1 | 0.02 | 118 | 0.02027 | 0.02027 |
| GO:0008614 | pyridoxine metabolic process | 3 | 1 | 0.02 | 119 | 0.02027 | 0.02027 |
| GO:0008615 | pyridoxine biosynthetic process | 3 | 1 | 0.02 | 120 | 0.02027 | 0.02027 |
| GO:0043567 | regulation of insulin-like growth factor… | 3 | 1 | 0.02 | 121 | 0.02027 | 0.02027 |
| GO:0045678 | positive regulation of R7 cell different… | 3 | 1 | 0.02 | 122 | 0.02027 | 0.02027 |